home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / multiroots / test_funcs.c < prev    next >
Encoding:
C/C++ Source or Header  |  2001-05-14  |  16.2 KB  |  745 lines
  1. /* multiroots/test_funcs.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <math.h>
  22. #include <gsl/gsl_vector.h>
  23. #include <gsl/gsl_matrix.h>
  24. #include <gsl/gsl_multiroots.h>
  25.  
  26. #include "test_funcs.h"
  27.  
  28. /* For information on testing see the following paper,
  29.  
  30.    J.J More, B.S. Garbow, K.E. Hillstrom, "Testing Unconstrained
  31.    Optimization Software", ACM Transactions on Mathematical Software,
  32.    Vol 7, No 1, (1981) p 17-41
  33.  
  34.    */
  35.  
  36. /* Rosenbrock Function */
  37.  
  38. gsl_multiroot_function_fdf rosenbrock =
  39. {&rosenbrock_f,
  40.  &rosenbrock_df,
  41.  &rosenbrock_fdf,
  42.  2, 0};
  43.  
  44. void
  45. rosenbrock_initpt (gsl_vector * x)
  46. {
  47.   gsl_vector_set (x, 0, -1.2);
  48.   gsl_vector_set (x, 1, 1.0);
  49. }
  50.  
  51. int
  52. rosenbrock_f (const gsl_vector * x, void *params, gsl_vector * f)
  53. {
  54.   double x0 = gsl_vector_get (x, 0);
  55.   double x1 = gsl_vector_get (x, 1);
  56.  
  57.   double y0 = 1 - x0;
  58.   double y1 = 10 * (x1 - x0 * x0);
  59.  
  60.   gsl_vector_set (f, 0, y0);
  61.   gsl_vector_set (f, 1, y1);
  62.  
  63.   params = 0;            /* avoid warning about unused parameters */
  64.  
  65.   return GSL_SUCCESS;
  66. }
  67.  
  68. int
  69. rosenbrock_df (const gsl_vector * x, void *params, gsl_matrix * df)
  70. {
  71.   double x0 = gsl_vector_get (x, 0);
  72.  
  73.   double df00 = -1;
  74.   double df01 = 0;
  75.   double df10 = -20 * x0;
  76.   double df11 = 10;
  77.  
  78.   gsl_matrix_set (df, 0, 0, df00);
  79.   gsl_matrix_set (df, 0, 1, df01);
  80.   gsl_matrix_set (df, 1, 0, df10);
  81.   gsl_matrix_set (df, 1, 1, df11);
  82.  
  83.   params = 0;            /* avoid warning about unused parameters */
  84.  
  85.   return GSL_SUCCESS;
  86. }
  87.  
  88. int
  89. rosenbrock_fdf (const gsl_vector * x, void *params,
  90.         gsl_vector * f, gsl_matrix * df)
  91. {
  92.   rosenbrock_f (x, params, f);
  93.   rosenbrock_df (x, params, df);
  94.  
  95.   return GSL_SUCCESS;
  96. }
  97.  
  98.  
  99. /* Freudenstein and Roth function */
  100.  
  101. gsl_multiroot_function_fdf roth =
  102. {&roth_f,
  103.  &roth_df,
  104.  &roth_fdf,
  105.  2, 0};
  106.  
  107. void
  108. roth_initpt (gsl_vector * x)
  109. {
  110.   gsl_vector_set (x, 0, 4.5);  /* changed from the value in the paper */
  111.   gsl_vector_set (x, 1, 3.5);  /* otherwise the problem is too hard */
  112. }
  113.  
  114. int
  115. roth_f (const gsl_vector * x, void *params, gsl_vector * f)
  116. {
  117.   double x0 = gsl_vector_get (x, 0);
  118.   double x1 = gsl_vector_get (x, 1);
  119.  
  120.   double y0 = -13.0 + x0 + ((5.0 - x1)*x1 - 2.0)*x1;
  121.   double y1 = -29.0 + x0 + ((x1 + 1.0)*x1 - 14.0)*x1;
  122.  
  123.   gsl_vector_set (f, 0, y0);
  124.   gsl_vector_set (f, 1, y1);
  125.  
  126.   params = 0;            /* avoid warning about unused parameters */
  127.  
  128.   return GSL_SUCCESS;
  129. }
  130.  
  131. int
  132. roth_df (const gsl_vector * x, void *params, gsl_matrix * df)
  133. {
  134.   double x1 = gsl_vector_get (x, 1);
  135.  
  136.   double df00 = 1;
  137.   double df01 = -3 * x1 * x1 + 10 * x1 - 2;
  138.   double df10 = 1;
  139.   double df11 = 3 * x1 * x1 + 2 * x1 - 14;
  140.  
  141.   gsl_matrix_set (df, 0, 0, df00);
  142.   gsl_matrix_set (df, 0, 1, df01);
  143.   gsl_matrix_set (df, 1, 0, df10);
  144.   gsl_matrix_set (df, 1, 1, df11);
  145.  
  146.   params = 0;            /* avoid warning about unused parameters */
  147.  
  148.   return GSL_SUCCESS;
  149. }
  150.  
  151. int
  152. roth_fdf (const gsl_vector * x, void *params,
  153.         gsl_vector * f, gsl_matrix * df)
  154. {
  155.   roth_f (x, params, f);
  156.   roth_df (x, params, df);
  157.  
  158.   return GSL_SUCCESS;
  159. }
  160.  
  161.  
  162.  
  163. /* Powell badly scaled function */
  164.  
  165. gsl_multiroot_function_fdf powellscal =
  166. {&powellscal_f,
  167.  &powellscal_df,
  168.  &powellscal_fdf,
  169.  2, 0};
  170.  
  171. void
  172. powellscal_initpt (gsl_vector * x)
  173. {
  174.   gsl_vector_set (x, 0, 0.0);
  175.   gsl_vector_set (x, 1, 1.0);
  176. }
  177.  
  178. int
  179. powellscal_f (const gsl_vector * x, void *params, gsl_vector * f)
  180. {
  181.   double x0 = gsl_vector_get (x, 0);
  182.   double x1 = gsl_vector_get (x, 1);
  183.  
  184.   double y0 = 10000.0 * x0 * x1 - 1.0;
  185.   double y1 = exp (-x0) + exp (-x1) - 1.0001;
  186.  
  187.   gsl_vector_set (f, 0, y0);
  188.   gsl_vector_set (f, 1, y1);
  189.  
  190.   params = 0;            /* avoid warning about unused parameters */
  191.  
  192.   return GSL_SUCCESS;
  193. }
  194.  
  195. int
  196. powellscal_df (const gsl_vector * x, void *params, gsl_matrix * df)
  197. {
  198.   double x0 = gsl_vector_get (x, 0);
  199.   double x1 = gsl_vector_get (x, 1);
  200.  
  201.   double df00 = 10000.0 * x1, df01 = 10000.0 * x0;
  202.   double df10 = -exp (-x0), df11 = -exp (-x1);
  203.  
  204.   gsl_matrix_set (df, 0, 0, df00);
  205.   gsl_matrix_set (df, 0, 1, df01);
  206.   gsl_matrix_set (df, 1, 0, df10);
  207.   gsl_matrix_set (df, 1, 1, df11);
  208.  
  209.   params = 0;            /* avoid warning about unused parameters */
  210.  
  211.   return GSL_SUCCESS;
  212. }
  213.  
  214. int
  215. powellscal_fdf (const gsl_vector * x, void *params,
  216.           gsl_vector * f, gsl_matrix * df)
  217. {
  218.   powellscal_f (x, params, f);
  219.   powellscal_df (x, params, df);
  220.  
  221.   return GSL_SUCCESS;
  222. }
  223.  
  224.  
  225. /* Brown badly scaled function */
  226.  
  227. gsl_multiroot_function_fdf brownscal =
  228. {&brownscal_f,
  229.  &brownscal_df,
  230.  &brownscal_fdf,
  231.  2, 0};
  232.  
  233. void
  234. brownscal_initpt (gsl_vector * x)
  235. {
  236.   gsl_vector_set (x, 0, 1.0);
  237.   gsl_vector_set (x, 1, 1.0);
  238. }
  239.  
  240. int
  241. brownscal_f (const gsl_vector * x, void *params, gsl_vector * f)
  242. {
  243.   double x0 = gsl_vector_get (x, 0);
  244.   double x1 = gsl_vector_get (x, 1);
  245.  
  246.   double y0 = x0 - 1e6;
  247.   double y1 = x0 * x1 - 2;
  248.  
  249.   gsl_vector_set (f, 0, y0);
  250.   gsl_vector_set (f, 1, y1);
  251.  
  252.   params = 0;            /* avoid warning about unused parameters */
  253.  
  254.   return GSL_SUCCESS;
  255. }
  256.  
  257. int
  258. brownscal_df (const gsl_vector * x, void *params, gsl_matrix * df)
  259. {
  260.   double x0 = gsl_vector_get (x, 0);
  261.   double x1 = gsl_vector_get (x, 1);
  262.  
  263.   double df00 = 1.0, df01 = 0.0;
  264.   double df10 = x1, df11 = x0;
  265.  
  266.   gsl_matrix_set (df, 0, 0, df00);
  267.   gsl_matrix_set (df, 0, 1, df01);
  268.   gsl_matrix_set (df, 1, 0, df10);
  269.   gsl_matrix_set (df, 1, 1, df11);
  270.  
  271.   params = 0;            /* avoid warning about unused parameters */
  272.  
  273.   return GSL_SUCCESS;
  274. }
  275.  
  276. int
  277. brownscal_fdf (const gsl_vector * x, void *params,
  278.           gsl_vector * f, gsl_matrix * df)
  279. {
  280.   brownscal_f (x, params, f);
  281.   brownscal_df (x, params, df);
  282.  
  283.   return GSL_SUCCESS;
  284. }
  285.  
  286.  
  287. /* Powell Singular Function */
  288.  
  289. gsl_multiroot_function_fdf powellsing =
  290. {&powellsing_f,
  291.  &powellsing_df,
  292.  &powellsing_fdf,
  293.  4, 0};
  294.  
  295. void
  296. powellsing_initpt (gsl_vector * x)
  297. {
  298.   gsl_vector_set (x, 0, 3.0);
  299.   gsl_vector_set (x, 1, -1.0);
  300.   gsl_vector_set (x, 2, 0.0);
  301.   gsl_vector_set (x, 3, 1.0);
  302. }
  303.  
  304. int
  305. powellsing_f (const gsl_vector * x, void *params, gsl_vector * f)
  306. {
  307.   double x0 = gsl_vector_get (x, 0);
  308.   double x1 = gsl_vector_get (x, 1);
  309.   double x2 = gsl_vector_get (x, 2);
  310.   double x3 = gsl_vector_get (x, 3);
  311.  
  312.   double y0 = x0 + 10 * x1;
  313.   double y1 = sqrt (5.0) * (x2 - x3);
  314.   double y2 = pow (x1 - 2 * x2, 2.0);
  315.   double y3 = sqrt (10.0) * pow (x0 - x3, 2.0);
  316.  
  317.   gsl_vector_set (f, 0, y0);
  318.   gsl_vector_set (f, 1, y1);
  319.   gsl_vector_set (f, 2, y2);
  320.   gsl_vector_set (f, 3, y3);
  321.  
  322.   params = 0;            /* avoid warning about unused parameters */
  323.  
  324.   return GSL_SUCCESS;
  325. }
  326.  
  327. int
  328. powellsing_df (const gsl_vector * x, void *params, gsl_matrix * df)
  329. {
  330.   double x0 = gsl_vector_get (x, 0);
  331.   double x1 = gsl_vector_get (x, 1);
  332.   double x2 = gsl_vector_get (x, 2);
  333.   double x3 = gsl_vector_get (x, 3);
  334.  
  335.   double df00 = 1, df01 = 10, df02 = 0, df03 = 0;
  336.   double df10 = 0, df11 = 0, df12 = sqrt (5.0), df13 = -df12;
  337.   double df20 = 0, df21 = 2 * (x1 - 2 * x2), df22 = -2 * df21, df23 = 0;
  338.   double df30 = 2 * sqrt (10.0) * (x0 - x3), df31 = 0, df32 = 0, df33 = -df30;
  339.  
  340.   gsl_matrix_set (df, 0, 0, df00);
  341.   gsl_matrix_set (df, 0, 1, df01);
  342.   gsl_matrix_set (df, 0, 2, df02);
  343.   gsl_matrix_set (df, 0, 3, df03);
  344.  
  345.   gsl_matrix_set (df, 1, 0, df10);
  346.   gsl_matrix_set (df, 1, 1, df11);
  347.   gsl_matrix_set (df, 1, 2, df12);
  348.   gsl_matrix_set (df, 1, 3, df13);
  349.  
  350.   gsl_matrix_set (df, 2, 0, df20);
  351.   gsl_matrix_set (df, 2, 1, df21);
  352.   gsl_matrix_set (df, 2, 2, df22);
  353.   gsl_matrix_set (df, 2, 3, df23);
  354.  
  355.   gsl_matrix_set (df, 3, 0, df30);
  356.   gsl_matrix_set (df, 3, 1, df31);
  357.   gsl_matrix_set (df, 3, 2, df32);
  358.   gsl_matrix_set (df, 3, 3, df33);
  359.  
  360.   params = 0;            /* avoid warning about unused parameters */
  361.  
  362.   return GSL_SUCCESS;
  363. }
  364.  
  365. int
  366. powellsing_fdf (const gsl_vector * x, void *params,
  367.             gsl_vector * f, gsl_matrix * df)
  368. {
  369.   powellsing_f (x, params, f);
  370.   powellsing_df (x, params, df);
  371.  
  372.   return GSL_SUCCESS;
  373. }
  374.  
  375.  
  376. /* Wood function */
  377.  
  378. gsl_multiroot_function_fdf wood =
  379. {&wood_f,
  380.  &wood_df,
  381.  &wood_fdf,
  382.  4, 0};
  383.  
  384. void
  385. wood_initpt (gsl_vector * x)
  386. {
  387.   gsl_vector_set (x, 0, -3.0);
  388.   gsl_vector_set (x, 1, -1.0);
  389.   gsl_vector_set (x, 2, -3.0);
  390.   gsl_vector_set (x, 3, -1.0);
  391. }
  392.  
  393. int
  394. wood_f (const gsl_vector * x, void *params, gsl_vector * f)
  395. {
  396.   double x0 = gsl_vector_get (x, 0);
  397.   double x1 = gsl_vector_get (x, 1);
  398.   double x2 = gsl_vector_get (x, 2);
  399.   double x3 = gsl_vector_get (x, 3);
  400.  
  401.   double t1 = x1 - x0 * x0;
  402.   double t2 = x3 - x2 * x2;
  403.  
  404.   double y0 = -200.0 * x0 * t1 - (1 - x0);
  405.   double y1 = 200.0 * t1 + 20.2 * (x1 - 1) + 19.8 * (x3 - 1);
  406.   double y2 = -180.0 * x2 * t2 - (1 - x2);
  407.   double y3 = 180.0 * t2 + 20.2 * (x3 - 1) + 19.8 * (x1 - 1);
  408.  
  409.   gsl_vector_set (f, 0, y0);
  410.   gsl_vector_set (f, 1, y1);
  411.   gsl_vector_set (f, 2, y2);
  412.   gsl_vector_set (f, 3, y3);
  413.  
  414.   params = 0;            /* avoid warning about unused parameters */
  415.  
  416.   return GSL_SUCCESS;
  417. }
  418.  
  419. int
  420. wood_df (const gsl_vector * x, void *params, gsl_matrix * df)
  421. {
  422.   double x0 = gsl_vector_get (x, 0);
  423.   double x1 = gsl_vector_get (x, 1);
  424.   double x2 = gsl_vector_get (x, 2);
  425.   double x3 = gsl_vector_get (x, 3);
  426.  
  427.   double t1 = x1 - 3 * x0 * x0;
  428.   double t2 = x3 - 3 * x2 * x2;
  429.  
  430.   double df00 = -200.0 * t1 + 1, df01 = -200.0 * x0, df02 = 0, df03 = 0;
  431.   double df10 = -400.0*x0, df11 = 200.0 + 20.2, df12 = 0, df13 = 19.8;
  432.   double df20 = 0, df21 = 0, df22 = -180.0 * t2 + 1, df23 = -180.0 * x2;
  433.   double df30 = 0, df31 = 19.8, df32 = -2 * 180 * x2, df33 = 180.0 + 20.2;
  434.  
  435.   gsl_matrix_set (df, 0, 0, df00);
  436.   gsl_matrix_set (df, 0, 1, df01);
  437.   gsl_matrix_set (df, 0, 2, df02);
  438.   gsl_matrix_set (df, 0, 3, df03);
  439.  
  440.   gsl_matrix_set (df, 1, 0, df10);
  441.   gsl_matrix_set (df, 1, 1, df11);
  442.   gsl_matrix_set (df, 1, 2, df12);
  443.   gsl_matrix_set (df, 1, 3, df13);
  444.  
  445.   gsl_matrix_set (df, 2, 0, df20);
  446.   gsl_matrix_set (df, 2, 1, df21);
  447.   gsl_matrix_set (df, 2, 2, df22);
  448.   gsl_matrix_set (df, 2, 3, df23);
  449.  
  450.   gsl_matrix_set (df, 3, 0, df30);
  451.   gsl_matrix_set (df, 3, 1, df31);
  452.   gsl_matrix_set (df, 3, 2, df32);
  453.   gsl_matrix_set (df, 3, 3, df33);
  454.  
  455.   params = 0;            /* avoid warning about unused parameters */
  456.  
  457.   return GSL_SUCCESS;
  458. }
  459.  
  460. int
  461. wood_fdf (const gsl_vector * x, void *params,
  462.             gsl_vector * f, gsl_matrix * df)
  463. {
  464.   wood_f (x, params, f);
  465.   wood_df (x, params, df);
  466.  
  467.   return GSL_SUCCESS;
  468. }
  469.  
  470.  
  471. /* Helical Valley Function */
  472.  
  473. gsl_multiroot_function_fdf helical =
  474. {&helical_f,
  475.  &helical_df,
  476.  &helical_fdf,
  477.  3, 0};
  478.  
  479. void
  480. helical_initpt (gsl_vector * x)
  481. {
  482.   gsl_vector_set (x, 0, -1.0);
  483.   gsl_vector_set (x, 1, 0.0);
  484.   gsl_vector_set (x, 2, 0.0);
  485. }
  486.  
  487. int
  488. helical_f (const gsl_vector * x, void *params, gsl_vector * f)
  489. {
  490.   double x0 = gsl_vector_get (x, 0);
  491.   double x1 = gsl_vector_get (x, 1);
  492.   double x2 = gsl_vector_get (x, 2);
  493.  
  494.   double t1, t2;
  495.   double y0, y1, y2;
  496.  
  497.   if (x0 > 0) 
  498.     {
  499.       t1 = atan(x1/x0) / (2.0 * M_PI);
  500.     }
  501.   else if (x0 < 0)
  502.     {
  503.       t1 = 0.5 + atan(x1/x0) / (2.0 * M_PI);
  504.     }
  505.   else
  506.     {
  507.       t1 = 0.25 * (x1 > 0 ? +1 : -1);
  508.     }
  509.  
  510.   t2 = sqrt(x0*x0 + x1*x1) ;
  511.   
  512.   y0 = 10 * (x2 - 10 * t1);
  513.   y1 = 10 * (t2 - 1);
  514.   y2 = x2 ;
  515.  
  516.   gsl_vector_set (f, 0, y0);
  517.   gsl_vector_set (f, 1, y1);
  518.   gsl_vector_set (f, 2, y2);
  519.  
  520.   params = 0;            /* avoid warning about unused parameters */
  521.  
  522.   return GSL_SUCCESS;
  523. }
  524.  
  525. int
  526. helical_df (const gsl_vector * x, void *params, gsl_matrix * df)
  527. {
  528.   double x0 = gsl_vector_get (x, 0);
  529.   double x1 = gsl_vector_get (x, 1);
  530.  
  531.   double t = x0 * x0 + x1 * x1 ;
  532.   double t1 = 2 * M_PI * t ;
  533.   double t2 = sqrt(t) ;
  534.  
  535.   double df00 = 100*x1/t1, df01 = -100.0 * x0/t1, df02 = 10.0;
  536.   double df10 = 10*x0/t2, df11 = 10*x1/t2, df12 = 0;
  537.   double df20 = 0, df21 = 0, df22 = 1.0;
  538.  
  539.   gsl_matrix_set (df, 0, 0, df00);
  540.   gsl_matrix_set (df, 0, 1, df01);
  541.   gsl_matrix_set (df, 0, 2, df02);
  542.  
  543.   gsl_matrix_set (df, 1, 0, df10);
  544.   gsl_matrix_set (df, 1, 1, df11);
  545.   gsl_matrix_set (df, 1, 2, df12);
  546.  
  547.   gsl_matrix_set (df, 2, 0, df20);
  548.   gsl_matrix_set (df, 2, 1, df21);
  549.   gsl_matrix_set (df, 2, 2, df22);
  550.  
  551.   params = 0;            /* avoid warning about unused parameters */
  552.  
  553.   return GSL_SUCCESS;
  554. }
  555.  
  556. int
  557. helical_fdf (const gsl_vector * x, void *params,
  558.             gsl_vector * f, gsl_matrix * df)
  559. {
  560.   helical_f (x, params, f);
  561.   helical_df (x, params, df);
  562.  
  563.   return GSL_SUCCESS;
  564. }
  565.  
  566.  
  567. /* Discrete Boundary Value Function */
  568.  
  569. #define N 10
  570.  
  571. gsl_multiroot_function_fdf dbv =
  572. {&dbv_f,
  573.  &dbv_df,
  574.  &dbv_fdf,
  575.  N, 0};
  576.  
  577. void
  578. dbv_initpt (gsl_vector * x)
  579. {
  580.   size_t i;
  581.   double h = 1.0 / (N + 1.0);
  582.  
  583.   for (i = 0; i < N; i++)
  584.     {
  585.       double t = (i + 1) * h;
  586.       double z = t * (t - 1);
  587.       gsl_vector_set (x, i, z);
  588.     }
  589. }
  590.  
  591. int
  592. dbv_f (const gsl_vector * x, void *params, gsl_vector * f)
  593. {
  594.   size_t i;
  595.  
  596.   double h = 1.0 / (N + 1.0);
  597.  
  598.   for (i = 0; i < N; i++)
  599.     {
  600.       double z, ti = (i + 1) * h;
  601.       double xi = 0, xim1 = 0, xip1 = 0;
  602.  
  603.       xi = gsl_vector_get (x, i);
  604.       
  605.       if (i > 0)
  606.         xim1 = gsl_vector_get (x, i - 1);
  607.  
  608.       if (i < N - 1)
  609.         xip1 = gsl_vector_get (x, i + 1);
  610.  
  611.       z = 2 * xi - xim1 - xip1 + h * h * pow(xi + ti + 1, 3.0) / 2.0;
  612.  
  613.       gsl_vector_set (f, i, z);
  614.  
  615.     }
  616.  
  617.   params = 0;            /* avoid warning about unused parameters */
  618.  
  619.   return GSL_SUCCESS;
  620. }
  621.  
  622. int
  623. dbv_df (const gsl_vector * x, void *params, gsl_matrix * df)
  624. {
  625.   size_t i, j;
  626.  
  627.   double h = 1.0 / (N + 1.0);
  628.  
  629.   for (i = 0; i < N; i++)
  630.     for (j = 0; j < N; j++)
  631.       gsl_matrix_set (df, i, j, 0.0);
  632.  
  633.   for (i = 0; i < N; i++)
  634.     {
  635.       double dz_dxi, ti = (i + 1) * h;
  636.  
  637.       double xi = gsl_vector_get (x, i);
  638.       
  639.       dz_dxi = 2.0 + (3.0 / 2.0) * h * h * pow(xi + ti + 1, 2.0) ;
  640.       
  641.       gsl_matrix_set (df, i, i, dz_dxi);
  642.  
  643.       if (i > 0)
  644.         gsl_matrix_set (df, i, i-1, -1.0);
  645.  
  646.       if (i < N - 1)
  647.         gsl_matrix_set (df, i, i+1, -1.0);
  648.  
  649.     }
  650.  
  651.   params = 0;            /* avoid warning about unused parameters */
  652.  
  653.   return GSL_SUCCESS;
  654. }
  655.  
  656. int
  657. dbv_fdf (const gsl_vector * x, void *params,
  658.             gsl_vector * f, gsl_matrix * df)
  659. {
  660.   dbv_f (x, params, f);
  661.   dbv_df (x, params, df);
  662.  
  663.   return GSL_SUCCESS;
  664. }
  665.  
  666. /* Trigonometric Function */
  667.  
  668. gsl_multiroot_function_fdf trig =
  669. {&trig_f,
  670.  &trig_df,
  671.  &trig_fdf,
  672.  N, 0};
  673.  
  674. void
  675. trig_initpt (gsl_vector * x)
  676. {
  677.   size_t i;
  678.  
  679.   for (i = 0; i < N; i++)       /* choose an initial point which converges */
  680.     {
  681.       gsl_vector_set (x, i, 0.05);   
  682.     }
  683. }
  684.  
  685. int
  686. trig_f (const gsl_vector * x, void *params, gsl_vector * f)
  687. {
  688.   size_t i;
  689.   double sum = 0;
  690.  
  691.   for (i = 0; i < N; i++)
  692.     {
  693.       sum += cos(gsl_vector_get(x,i));
  694.     }
  695.  
  696.   for (i = 0; i < N; i++)
  697.     {
  698.       double xi = gsl_vector_get (x,i);
  699.       double z = N - sum + (i + 1) * (1 - cos(xi)) - sin(xi);
  700.  
  701.       gsl_vector_set (f, i, z);
  702.     }
  703.  
  704.   params = 0;            /* avoid warning about unused parameters */
  705.  
  706.   return GSL_SUCCESS;
  707. }
  708.  
  709. int
  710. trig_df (const gsl_vector * x, void *params, gsl_matrix * df)
  711. {
  712.   size_t i, j;
  713.  
  714.   for (i = 0; i < N; i++)
  715.     {
  716.       for (j = 0; j < N; j++)
  717.         {
  718.           double dz;
  719.           double xi = gsl_vector_get(x, i);
  720.           double xj = gsl_vector_get(x, j);
  721.  
  722.           if (j == i)
  723.             dz = sin(xi) + (i + 1) * sin(xi) - cos(xi);
  724.           else
  725.             dz = sin(xj);
  726.           
  727.           gsl_matrix_set(df, i, j, dz);
  728.         }
  729.     }
  730.  
  731.   params = 0;            /* avoid warning about unused parameters */
  732.  
  733.   return GSL_SUCCESS;
  734. }
  735.  
  736. int
  737. trig_fdf (const gsl_vector * x, void *params,
  738.             gsl_vector * f, gsl_matrix * df)
  739. {
  740.   trig_f (x, params, f);
  741.   trig_df (x, params, df);
  742.  
  743.   return GSL_SUCCESS;
  744. }
  745.